#include <opencv2/opencv.hpp>
#include <list>
#include "tools.h"


#define PI 3.1416

using namespace cv;
using namespace std;

static const char _FILE_PREFIX[] = "out/img";
static const long long primeNum[] = {
        104743, 105023, 105359, 105613, 104759, 105031, 105361, 105619, 104761,
        105037, 105367, 105649, 104773, 105071, 105373, 105653, 104779, 105097,
        105379, 105667, 104789, 105107, 105389, 105673, 104801, 105137, 105397,
        105683, 104803, 105143, 105401, 105691, 104827, 105167, 105407, 105701,
        104831, 105173, 105437, 105727, 104849, 105199, 105449, 105733, 104851,
        105211, 105467, 105751, 104869, 105227, 105491, 105761, 104879, 105229,
        105499, 105767, 104891, 105239, 105503, 105769, 104911, 105251, 105509,
        105817, 104917, 105253, 105517, 105829, 104933, 105263, 105527, 105863,
        104947, 105269, 105529, 105871, 104953, 105277, 105533, 105883, 104959,
        105319, 105541, 105899, 104971, 105323, 105557, 105907, 104987, 105331,
        105563, 105913, 104999, 105337, 105601, 105929, 105019, 105341, 105607,
        105943
};

/*///////////////////////////////////////////////////////////////////////////////
 *  Vessel extract related functions
 */
Mat find_areaRange(Mat img, int minarea, int maxarea);

Mat find_connected(Mat img, int area = 1000);

Mat getLineMask(int theta, int masksize);

Mat rotatex(Mat mask);

Mat getBaseMask(int theta, int masksize);

Mat drawLine(cv::Point pt0, cv::Point pt1, Mat orgmask);

Mat getLinePts(cv::Point pt0, cv::Point pt1);

Mat grow_boundary(Mat img, Mat mask, int erosionsize = 21, int iterations = 20);

Mat global_normalize(Mat img, Mat mask);

Mat get_lineresponse(Mat img, int linescale, int w = 15);

Mat create_mask(Mat _img, double thresh = 0.001);

Mat vessel_seg(Mat img);

/*
 *  This function returen a mask indicating the optic disk and the center
 *  bright area that should not be regarded as area of leakage.
 *
 *  Retuned mask is a 1-chnnel and cv_8u type image.
 */
Mat getOpticDisk(Mat input);

/*
 *  This function performs local normalization, which is good for now.
 */
Mat localNorm(Mat image, float sigma1, float sigma2, InputArray mask);

Mat find_roundObject(Mat img, int radius);

/*
 *  This function return CV_8UC1 Data
 */
Mat find_connected(Mat img, int area) {
    Mat label, stats, centroids;
    connectedComponentsWithStats(img, label, stats, centroids, 4, CV_32S);
    Mat oimg(img.size[0], img.size[1], CV_8UC1);
    oimg = 0;
    for (int i = 1; i < stats.size[0]; i++) {
        if (stats.at<int32_t>(i, CC_STAT_AREA) < area)
            continue;
        for (int x = stats.at<int32_t>(i, CC_STAT_TOP);
             x < stats.at<int32_t>(i, CC_STAT_TOP) + stats.at<int32_t>(i, CC_STAT_HEIGHT);
             x++) {
            for (int y = stats.at<int32_t>(i, CC_STAT_LEFT);
                 y < stats.at<int32_t>(i, CC_STAT_LEFT) + stats.at<int32_t>(i, CC_STAT_WIDTH);
                 y++) {
                if (label.at<uint32_t>(x, y) == i) {
                    oimg.at<uchar>(x, y) = 255;
                }
            }
        }
//        printf("handled %d, area %d\n", i, stats.at<int32_t>(i,CC_STAT_AREA));
    }
    return oimg;
}

Mat find_areaRange(Mat img, int minarea, int maxarea) {
    Mat label, stats, centroids;
    connectedComponentsWithStats(img, label, stats, centroids, 8, CV_32S);
    Mat oimg(img.size[0], img.size[1], CV_8UC1);
    oimg = 0;
    for (int i = 1; i < stats.size[0]; i++) {
        if (stats.at<int32_t>(i, CC_STAT_AREA) < minarea ||
            stats.at<int32_t>(i, CC_STAT_AREA) > maxarea)
            continue;
        for (int x = stats.at<int32_t>(i, CC_STAT_TOP);
             x < stats.at<int32_t>(i, CC_STAT_TOP) + stats.at<int32_t>(i, CC_STAT_HEIGHT);
             x++) {
            for (int y = stats.at<int32_t>(i, CC_STAT_LEFT);
                 y < stats.at<int32_t>(i, CC_STAT_LEFT) + stats.at<int32_t>(i, CC_STAT_WIDTH);
                 y++) {
                if (label.at<uint32_t>(x, y) == i) {
                    oimg.at<uchar>(x, y) = 255;
                }
            }
        }
//        printf("handled %d, area %d\n", i, stats.at<int32_t>(i,CC_STAT_AREA));
    }
    return oimg;
}

Mat vessel_seg(Mat img) {
    Mat mask_img = create_mask(img);
    vector<Mat> image_channels;
    split(img, image_channels);
    Mat gChannel = image_channels[1];
//    medianBlur(gChannel,gChannel,3);
//    Mat ignore = FindReflection(image_channels[0], mask_img);
    // split(mask_img, mask_channels);
    // Mat mask_gChannel = mask_channels[1];
    Mat mask_gChannel = mask_img;
    gChannel = 255 - gChannel;    /* Í¼Ïñ»Ò¶È·´×ª£¬Ê¹Ñª¹Ü»Ò¶È¸ßÓÚ±³¾° */
    Mat dilated = grow_boundary(gChannel, mask_gChannel);
    Mat dilated_double = Mat::zeros(mask_img.rows, mask_img.cols, CV_64FC1);
    Mat dilated_norm = Mat::zeros(mask_img.rows, mask_img.cols, CV_64FC1);
    dilated.convertTo(dilated_double, CV_64FC1, 1.0 / 255.0);
    dilated_norm = global_normalize(dilated_double, mask_gChannel);
    Mat line_response = Mat::zeros(mask_img.rows, mask_img.cols, CV_64FC1);
    Mat line_response_norm = Mat::zeros(mask_img.rows, mask_img.cols, CV_64FC1);

    int scales = 0;
    int w = 21;
    int step = 2;

#pragma omp parallel for shared(dilated_norm) private(line_response,line_response_norm)
    for (int line = 1; line < 17; line += step) {
        line_response = get_lineresponse(dilated_double, line, w);
        line_response_norm = global_normalize(line_response, mask_gChannel);
        dilated_norm = dilated_norm + line_response_norm;
        scales++;
    }
    Mat avg_response = dilated_norm / (1 + scales);

    Mat vessels = Mat::zeros(img.rows, img.cols, CV_8UC1);
    double thresh = 0.48;
    for (int i = 0; i < img.rows; i++)
        for (int j = 0; j < img.cols; j++)
            if (avg_response.at<double>(i, j) > thresh) {
                vessels.at<uchar>(i, j) = 255;
            }
//    ignore = 255 - ignore;
//    erode(ignore, ignore, cv::Mat::ones(21, 21, CV_8U));
//    vessels = vessels.mul(ignore / 255);
    vessels = find_connected(vessels, 200);
    return vessels;
}

Mat get_lineresponse(Mat img, int linescale, int w) {
    Mat avg_kernel(w, w, CV_64FC1, Scalar(1));
    avg_kernel = avg_kernel / (w * w);
    Mat avg_response;
    Mat linemask(linescale, linescale, CV_8UC1, Scalar(0));
    Mat linemask_double(linescale, linescale, CV_64FC1, Scalar(0));
    Mat line_strength;
    Mat line_response = (-200.0) * Mat::ones(img.rows, img.cols, CV_64FC1);
    filter2D(img, avg_response, img.depth(), avg_kernel);
    for (int theta = 0; theta < 180; theta += 15) {
        linemask = getLineMask(theta, linescale) * 1.0;
        linemask.convertTo(linemask_double, CV_64FC1, 1.0 / 255.0);
        linemask_double = linemask_double / countNonZero(linemask_double);
        filter2D(img, line_strength, img.depth(), linemask_double);
        line_strength = line_strength - avg_response;
        for (int i = 0; i < img.rows; i++)
            for (int j = 0; j < img.cols; j++) {
                if (line_strength.at<double>(i, j) > line_response.at<double>(i, j))
                    line_response.at<double>(i, j) = line_strength.at<double>(i, j);
            }
    }
    return line_response;
}

Mat grow_boundary(Mat img, Mat mask_img, int erosionsize, int iterations) {
    Mat dilated = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
    for (int i = 0; i < mask_img.rows; i++)
        for (int j = 0; j < mask_img.cols; j++) {
            if (i == 0 || i == mask_img.rows - 1)
                mask_img.at<uchar>(i, j) = 0;
            if (j == 0 || j == mask_img.cols - 1)
                mask_img.at<uchar>(i, j) = 0;
        }
    Mat erode_element = getStructuringElement(MORPH_ELLIPSE, Size(erosionsize, erosionsize));
    Mat erode_mask = mask_img.clone();
    erode(mask_img, erode_mask, erode_element);
    dilated = img.mul(erode_mask / 255);
    Mat oldmask = erode_mask.clone();
    int filter_rows[9] = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
    int filter_cols[9] = {-1, -1, -1, 0, 0, 0, 1, 1, 1};
    Mat newmask = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
    Mat outerborder = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
    int pixelrow, pixelcol;
    int total, count;

    Mat diamond(3, 3, CV_8U, Scalar(1));
    diamond.at<uchar>(0, 0) = 0;
    diamond.at<uchar>(2, 0) = 0;
    diamond.at<uchar>(0, 2) = 0;
    diamond.at<uchar>(2, 2) = 0;

    for (int i = 0; i < iterations; i++) {
        outerborder = Mat::zeros(mask_img.rows, mask_img.cols, CV_8UC1);
        dilate(oldmask, newmask, diamond);
        outerborder = newmask - oldmask;
        for (int j = 0; j < outerborder.rows; j++)
            for (int k = 0; k < outerborder.cols; k++)
                if (outerborder.at<uchar>(j, k)) {
                    total = 0;
                    count = 0;
                    for (int h = 0; h < 9; h++) {
                        pixelrow = j + filter_rows[h];
                        pixelcol = k + filter_cols[h];
                        if (pixelrow < outerborder.rows && pixelrow >= 0 \
 && pixelcol < outerborder.cols && pixelcol >= 0 \
 && oldmask.at<uchar>(pixelrow, pixelcol)) {
                            total = total + dilated.at<uchar>(pixelrow, pixelcol);
                            count++;
                        }
                    }
                    dilated.at<uchar>(j, k) = (uchar) (total / count);
                }
        oldmask = newmask.clone();
    }
    return dilated;
}

Mat global_normalize(Mat img, Mat mask) {
    Mat norm_img(img.rows, img.cols, CV_64FC1, Scalar(0));
    Scalar mean, stdDev;
    int used_pixels = countNonZero(mask);
    double diff;
    meanStdDev(img, mean, stdDev, mask);
    stdDev.val[0] = stdDev.val[0] * std::sqrt(used_pixels / (used_pixels - 1.0));
    for (int i = 0; i < norm_img.rows; i++)
        for (int j = 0; j < norm_img.cols; j++)
            if (mask.at<uchar>(i, j)) {
                diff = (img.at<double>(i, j) - mean.val[0]);
                if (diff > 0) {
                    norm_img.at<double>(i, j) = diff / stdDev.val[0];
                }

            }
    return norm_img;
}


Mat getLinePts(Point pt0, Point pt1) {
    int x0 = 0, x1 = 0, y0 = 0, y1 = 0;
    if (pt0.x < pt1.x) {
        x0 = pt0.x;
        y0 = pt0.y;
        x1 = pt1.x;
        y1 = pt1.y;
    } else {
        x0 = pt1.x;
        y0 = pt1.y;
        x1 = pt0.x;
        y1 = pt0.y;
    }
    int dx = x1 - x0;
    int dy = y1 - y0;
    int x = 0, y = 0;
    int maxi = 0, mini = 0, maxc = 0;
    maxc = (abs(dx) > abs(dy)) ? abs(dx) : abs(dy);
    int ind = 0;//ÒòÎªÏÂ±ê´Ó0¿ªÊ¼£»
    uchar *dataPtr;
    Mat lineMask(maxc + 1, 2, CV_8U, Scalar(0));
    if (dx == 0) {
        x = x0;
        //if (dy < 0)
        {
            maxi = (y0 > y1) ? y0 : y1;
            mini = (y0 < y1) ? y0 : y1;
        }
        for (y = mini; y <= maxi; y++) {
            dataPtr = lineMask.ptr<uchar>(ind);
            dataPtr[0] = (uchar) y;
            dataPtr[1] = (uchar) x;
            ind++;
        }
    } else if (abs(dy) > abs(dx)) {
        //if (dy < 0)
        {
            maxi = (y0 > y1) ? y0 : y1;
            mini = (y0 < y1) ? y0 : y1;
        }
        for (y = mini; y <= maxi; y++) {
            x = (int) (((dx + 0.0000) / dy + 0.0000) * (y - y0) + x0);
            dataPtr = lineMask.ptr<uchar>(ind);
            dataPtr[0] = (uchar) y;
            dataPtr[1] = (uchar) x;
            ind++;
        }
    } else {
        for (x = x0; x <= x1; x++) {
            y = (int) (((dy + 0.0000) / dx + 0.0000) * (x - x0) + y0);
            dataPtr = lineMask.ptr<uchar>(ind);
            dataPtr[0] = (uchar) y;
            dataPtr[1] = (uchar) x;
            ind++;
        }
    }
    return lineMask;
}

Mat drawLine(Point pt0, Point pt1, Mat orgmask) {
    Mat mask = orgmask.clone();
    int maxc = 0;
    int dx = pt0.x - pt1.x;
    int dy = pt0.y - pt1.y;
    maxc = (abs(dx) > abs(dy)) ? abs(dx) : abs(dy);
    Mat tmpMat(maxc + 1, 2, CV_8U, Scalar(0));
    int tmp = 0, tmp1 = 0;
    tmpMat = getLinePts(pt0, pt1);
    uchar *dataPtr;
    uchar *dataPtr1;
    for (int i = 0; i < maxc + 1; i++) {
        dataPtr1 = tmpMat.ptr<uchar>(i);
        tmp = dataPtr1[0];
        dataPtr = mask.ptr<uchar>(tmp);
        tmp1 = dataPtr1[1];
        dataPtr[tmp1] = 255;
    }

    return mask;
}


Mat getBaseMask(int theta, int masksize) {
    Mat mask(masksize, masksize, CV_8U, Scalar(0));
    int halfsize = (masksize - 1) / 2;
    int x0 = 0, y0 = 0, x1 = 0, y1 = 0;
    Point pt0, pt1;
    uchar *dataPtr;
    if (theta == 0) {
        dataPtr = mask.ptr<uchar>(halfsize);
        for (int i = 0; i < mask.cols; i++)
            dataPtr[i] = 255;
    } else if (theta == 90) {
        for (int i = 0; i < mask.rows; i++) {
            dataPtr = mask.ptr<uchar>(i);
            dataPtr[halfsize] = 255;
        }
    } else {
        x0 = -halfsize;
        y0 = (int) (x0 * (tan(((float) theta + 0.0000) / 180 * PI) + 0.0000));
        if (y0 < -halfsize) {
            y0 = -halfsize;
            x0 = (int) (y0 * ((1 / (tan((float) theta * PI / 180) + 0.0000) + 0.0000)));
        }
        x1 = halfsize;
        y1 = (int) (x1 * (tan((float) theta * PI / 180) + 0.0000));
        if (y1 > halfsize) {
            y1 = halfsize;
            x1 = (int) (y1 * ((1 / (tan((float) theta * PI / 180) + 0.0000) + 0.0000)));
        }

        pt0.y = halfsize - y0;
        pt0.x = halfsize + x0;
        pt1.y = halfsize - y1;
        pt1.x = halfsize + x1;
        mask = drawLine(pt0, pt1, mask);
    }
    return mask;
}

Mat rotatex(Mat mask) {
    int h = 0, w = 0;
    h = mask.rows;
    w = mask.cols;
    Mat rotatedMask(h, w, CV_8U, Scalar(0));
    uchar *dataPtr, *dataPtr2;
    for (int i = 0; i < h; i++) {
        for (int j = 0; j < w; j++) {
            dataPtr = mask.ptr<uchar>(i);
            dataPtr2 = rotatedMask.ptr<uchar>(i);
            dataPtr2[j] = dataPtr[w - j - 1];
        }
    }
    return rotatedMask;
}


Mat getLineMask(int theta, int masksize) {
    Mat mask;
    Mat linemask;
    if (theta > 90) {
        mask = getBaseMask(180 - theta, masksize);
        linemask = rotatex(mask);
    } else {
        linemask = getBaseMask(theta, masksize);
    }
    return linemask;
}

/*
 * Imput image should be a 3-channel image
 */
Mat create_mask(Mat _img, double thresh) {
    CHECK(_img.channels() == 3);
    Mat img = myIm2Double(_img);
    vector<Mat> image_channels;
    split(img, image_channels);
    Mat rChannel(_img.rows, _img.cols, CV_32FC1, Scalar(0));
    for (int i = 0; i < 3; i++) {
        rChannel += image_channels[i].mul(image_channels[i]);
    }
    double maxValue, minValue;
    minMaxLoc(rChannel, &minValue, &maxValue);
    Mat seed = rChannel < thresh * maxValue;
    seed = 255 - seed;
    Mat blur;
    medianBlur(seed, blur, 11);
    return blur;
}

///////////////////////////////////////////////////////////////////


string getTimeString(void) {
    time_t rawtime;
    time(&rawtime);
    tm *time_ptr = localtime(&rawtime);
    char time_str[100];
    unsigned long len = strftime(time_str, sizeof(time_str), "%Y-%m-%d", time_ptr);
    return string(time_str, len);
}

static string _init_dir_path() {
    string filepath = string(_FILE_PREFIX) + "/" + getTimeString();
    boost::filesystem::path dir(filepath);
    if (!boost::filesystem::create_directories(dir)) {
        LOG(ERROR) << "cannot create dir : " << filepath;
    }
    return filepath;
}

string dir_path = _init_dir_path();

void genImgName(const char *name, char *oName, int index = -1) {
    if (index < 0)
        sprintf(oName, "%s/%s.png", dir_path.c_str(), name);
    else
        sprintf(oName, "%s/Fig%04d.%s.png", dir_path.c_str(), index, name);
}

void myImShow(Mat img, const char *name) {
    static int counter = 0;
    static char buf[100];
    if (name == NULL) {
        name = buf;
        sprintf(buf, "Figure %d", counter++);
    }
    namedWindow(name, CV_WINDOW_NORMAL);
    imshow(name, img);
    resizeWindow(name, 800, 600);
}

int myPause(char exitKey) {
    while (1) {
        char key = (char) cvWaitKey(0);
        if (key == -1)
            return key;
        else if (key == 27) // esc
            return key;
        else if (key == exitKey)
            return key;
    }
}


void myImWrite(Mat img, const char *name, int index) {
    int counter = 0;
    char buf[100];
    char numBuf[20];
    if (name == NULL) {
        sprintf(numBuf, "Figure_%d", counter);
        genImgName(numBuf, buf, index);
    } else {
        genImgName(name, buf, index);
    }
    imwrite(buf, myIm2Int(img));
}


double max(const Mat &m) {
    double vMax, vMin;
    minMaxLoc(m, &vMax, &vMin);
    return vMax;
}

double min(const Mat &m) {
    double vMax, vMin;
    minMaxLoc(m, &vMin, &vMax);
    return vMin;
}

Mat myImRead(const char *name, int flags, int index) {
    char nameBuf[100];
//    if(!myImExist(name)){
//        printf("Error: image not exist\n");
//        return NULL;
//    }
    genImgName(name, nameBuf, index);
    return imread(nameBuf, flags);
}


bool myImExist(const char *name) {
    char nameBuf[100];
    genImgName(name, nameBuf);
    return myFileExist(nameBuf);
}


bool myFileExist(const char *name) {
    struct stat buffer;
    return (stat(name, &buffer) == 0);
}

void myImDel(const char *name) {
    char nameBuf[100];
    if (!myImExist(name)) {
        printf("Warning: file not exist\n");
    }
    genImgName(name, nameBuf);
    remove(nameBuf);
}


Mat myIm2Double(Mat in) {
    Mat result;
    if (in.depth() == CV_32F) {
        return in;
    }
    CHECK_EQ(in.depth(), CV_8U);
    in.convertTo(result, CV_32FC(in.channels()), 1.0 / 255.0);
    return result;
}


Mat myIm2Int(Mat in) {
    Mat result;
    if (in.depth() == CV_8U) {
        return in;
    }
    CHECK(in.depth() == CV_32F || in.depth() == CV_64F);
    in.convertTo(result, CV_8UC(in.channels()), 255.0);
    return result;
}


/////////////////////////////////////////////////////////leakage
Mat elimate(Mat img, Mat mask) {
    CHECK(mask.depth() == CV_8U);
    CHECK(img.depth() == CV_8U);
    CHECK(mask.channels() == 1);
    CHECK(img.channels() == 1);
    Mat label, stats, centroids;
    connectedComponentsWithStats(img, label, stats, centroids, 8, CV_32S);
    Mat oimg(img.size[0], img.size[1], CV_8UC1);
    oimg = 0;
    int allarea = 0, area;
    for (int i = 1; i < stats.size[0]; i++) {
        allarea += stats.at<int32_t>(i, CC_STAT_AREA);
    }
    CHECK(stats.size[0] - 1 > 0) << "no points at all" << Common::printTrace();
    area = allarea / (stats.size[0] - 1);
    for (int i = 1; i < stats.size[0]; i++) {
        if (stats.at<int32_t>(i, CC_STAT_AREA) > area * 5)
            continue;
        int xmin = stats.at<int32_t>(i, CC_STAT_TOP);
        int xmax = xmin + stats.at<int32_t>(i, CC_STAT_HEIGHT);
        int ymin = stats.at<int32_t>(i, CC_STAT_LEFT);
        int ymax = ymin + stats.at<int32_t>(i, CC_STAT_WIDTH);
        bool breakout = false;
        for (int x = xmin; x < xmax; x++) {
            for (int y = ymin; y < ymax; y++) {
                if (mask.at<uchar>(x, y) != 0) {
                    breakout = true;
                }
            }
        }
        if (breakout)
            continue;
        for (int x = xmin; x < xmax; x++) {
            for (int y = ymin; y < ymax; y++) {
                if (label.at<uint32_t>(x, y) == i) {
                    oimg.at<uchar>(x, y) = 255;
                }
            }
        }
    }
    return oimg;
}


Mat FindReflection(Mat Rchannel, Mat mask) {
    CHECK(mask.depth() == CV_8U);
    CHECK(Rchannel.depth() == CV_8U);
    CHECK(Rchannel.channels() == 1);
    CHECK(mask.channels() == 1);
    Mat oimg;
    adaptiveThreshold(Rchannel, oimg, 255, ADAPTIVE_THRESH_MEAN_C, THRESH_BINARY, 151, -5);
    medianBlur(oimg, oimg, 7);
//    Mat mask1 = mask / 255;
//    Mat rr = myIm2Double(Rchannel);
//    double average = sum(rr.mul(myIm2Double(mask))).val[0] / sum(mask1).val[0];
//    return Rchannel > average * 255 * 1.5;
    return oimg;
}

Mat leakage_preprocess(Mat input) {
//    myImWrite(input, "oriimg", id);
    Mat mask = create_mask(input);
    Mat resizedImg;
    vector<Mat> channels;
    split(input, channels);
    Mat opticDisk = getOpticDisk(input);
    Mat nerveArea = FindReflection(channels[0], mask);
    resizedImg = grow_boundary(channels[1], mask, 60, 180);
    Mat normImg = localNorm(myIm2Double(resizedImg), 60, 300, mask);
    normImg = myIm2Int(normImg);
//    normalize(normImg,normImg,0,255,NORM_MINMAX, CV_8U,mask);
    Mat imgBG;
    medianBlur(normImg, imgBG, 101);
    Mat nobg = myIm2Double(normImg) - myIm2Double(imgBG);
    nobg *= 255;
    Mat nobg2;

    nobg = nobg > 25;
//    cv::AdaptiveThresholdTypes;
//    adaptiveThreshold(nobg,nobg,
//                      255,ADAPTIVE_THRESH_GAUSSIAN_C,
//                      THRESH_BINARY,201,-8);
    Mat kernel = getStructuringElement(MORPH_CROSS, Size(9, 9));
//    morphologyEx(nobg, nobg2, MORPH_OPEN, kernel);

    medianBlur(nobg, nobg2, 7);
//    nobg2 = nobg;

    Mat closeOp, openOp;
    morphologyEx(nobg2, closeOp, MORPH_CLOSE, kernel);
    morphologyEx(closeOp, openOp, MORPH_OPEN, kernel);
    openOp = closeOp;
    nobg2 = openOp;
    Mat elimateMask = opticDisk | nerveArea | (255 - mask);
    // nobg2 = elimate(nobg2, elimateMask);
    return nobg2;
}

void drawHist(Mat image, InputArray mask = cv::noArray(), int bins = 60) {
    const int channels[] = {0};
    Mat hist;
    float lranges[] = {0, 1};
    float *ranges[] = {lranges};
    const float **rr = (const float **) ranges;
    const int histSize[] = {bins};
    calcHist(&image, 1, channels, mask, // do not use mask
             hist, 1, histSize, rr,
             true, // the histogram is uniform
             false);

    int sbins = bins;
    int hbins = 20;
    int scale = 10;
    Mat histImg = Mat::zeros(hbins * 10, sbins * scale, CV_8UC3);
    double maxVal = 0;
    minMaxLoc(hist, 0, &maxVal, 0, 0);

    for (int s = 0; s < sbins; s++) {
        float binVal = hist.at<float>(s);
        int intensity = 255;
        rectangle(histImg, Point(s * scale, hbins * scale - 1),
                  Point((s + 1) * scale - 1, (int) (hbins * scale - hbins * binVal + 1)),
                  Scalar::all(intensity),
                  CV_FILLED);
    }

    namedWindow("H-S Histogram", 1);
    imshow("H-S Histogram", histImg);
    waitKey();
}


asd bleed_process(Mat input) {
//    myImWrite(input, "oriimg", id);
    vector<Mat> channels;
    split(input, channels);
    Mat mask = (create_mask(input));
    Mat nerveArea = FindReflection(channels[0], mask);
    Mat opticDisk = getOpticDisk(input);

    mask = myIm2Double(mask);
    Mat vessel = vessel_seg(input);
    dilate(vessel, vessel, Mat());
    dilate(vessel, vessel, Mat());

    dilate(vessel, vessel, Mat());
    Mat resizedImg;
//    split(input, channels);
//    Mat opticDisk = getOpticDisk(input);
//    Mat nerveArea = FindReflection(channels[0], mask);
    Mat tmp;
    input.convertTo(tmp, CV_32FC(input.channels()));
//    Mat c1 = myIm2Double(channels[1]);
    Mat c2;
    cv::cvtColor(tmp, c2, CV_BGR2HSV);
    split(c2, channels);
    Mat c1 = channels[0];
    c1 = c1.mul(mask);
//    drawHist(c1,myIm2Int(mask));
//    c1 = c1.mul(mask);
    double size = sum(mask)[0];
    double ave = sum(c1)[0] / size;
    LOG(INFO) << ave;
    // // imshow("1",c1);
    // // imshow("2",c1 < ave * 0.5);
    // // cvWaitKey(0);
    c2 = c1 < ave * 0.7;
    c2 = c2.mul(255 - vessel);
    Mat c3 = c1 > ave * 1.4;
    Mat elimateMask = opticDisk | nerveArea | (255 - myIm2Int(mask));
    c3 = elimate(c3, elimateMask);

    asd retval = {
            .bleed = c2.mul(myIm2Int(mask)),
            .exclude = c3.mul(myIm2Int(mask))
    };
    return retval;
//    resizedImg = grow_boundary(channels[1], mask, 60, 180);
//    Mat normImg = localNorm(myIm2Double(resizedImg), 60, 300, mask);
//    normImg = myIm2Int(normImg);
////    normalize(normImg,normImg,0,255,NORM_MINMAX, CV_8U,mask);
//    Mat imgBG;
//    medianBlur(normImg, imgBG, 31);
//    Mat nobg = myIm2Double(normImg) - myIm2Double(imgBG);
//
//    medianBlur(nobg, nobg, 3);
//    medianBlur(nobg, nobg, 3);
//    nobg = nobg < -15.0 / 255.0;
//
//    nobg = nobg.mul(255 - vessel);
//
//    medianBlur(nobg, nobg, 3);
//
//    return nobg;
}

double NormInv(double probability) {
    const double a1 = -39.6968302866538;
    const double a2 = 220.946098424521;
    const double a3 = -275.928510446969;
    const double a4 = 138.357751867269;
    const double a5 = -30.6647980661472;
    const double a6 = 2.50662827745924;

    const double b1 = -54.4760987982241;
    const double b2 = 161.585836858041;
    const double b3 = -155.698979859887;
    const double b4 = 66.8013118877197;
    const double b5 = -13.2806815528857;

    const double c1 = -7.78489400243029E-03;
    const double c2 = -0.322396458041136;
    const double c3 = -2.40075827716184;
    const double c4 = -2.54973253934373;
    const double c5 = 4.37466414146497;
    const double c6 = 2.93816398269878;

    const double d1 = 7.78469570904146E-03;
    const double d2 = 0.32246712907004;
    const double d3 = 2.445134137143;
    const double d4 = 3.75440866190742;

    //Define break-points
    // using Epsilon is wrong; see link above for reference to 0.02425 value
    //const double pLow = double.Epsilon;
    const double pLow = 0.02425;

    const double pHigh = 1 - pLow;

    //Define work variables
    double q;
    double result = 0;

    // if argument out of bounds.
    // set it to a value within desired precision.
    if (probability <= 0)
        probability = pLow;

    if (probability >= 1)
        probability = pHigh;

    if (probability < pLow) {
        //Rational approximation for lower region
        q = sqrt(-2 * log(probability));
        result = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) /
                 ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
    } else if (probability <= pHigh) {
        //Rational approximation for lower region
        q = probability - 0.5;
        double r = q * q;
        result = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q /
                 (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1);
    } else if (probability < 1) {
        //Rational approximation for upper region
        q = sqrt(-2 * log(1 - probability));
        result = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) /
                 ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
    }

    return result;
}

double NormInv(double probability, double mean, double sigma) {
    double x = NormInv(probability);
    return sigma * x + mean;
}


/*
 *  This function performs local normalization, which is good for now.
 */
Mat localNorm(Mat image, float sigma1, float sigma2, InputArray mask) {
    Mat floatGray, blurred1, blurred2, temp1, temp2, res;

//    image.convertTo(floatGray, CV_32FC1);
    CHECK(image.channels() == 1);
    CHECK(image.depth() == CV_32F);
    image.copyTo(floatGray);

    int blur1 = (int) (2 * ceil(-NormInv(0.05, 0, sigma1)) + 1);
    cv::GaussianBlur(floatGray, blurred1, cv::Size(blur1, blur1), sigma1);
    temp1 = floatGray - blurred1;

    cv::pow(temp1, 2.0, temp2);
    int blur2 = (int) (2 * ceil(-NormInv(0.05, 0, sigma2)) + 1);
    cv::GaussianBlur(temp2, blurred2, cv::Size(blur2, blur2), sigma2);
    cv::pow(blurred2, 0.5, temp2);

    floatGray = temp1 / temp2;
//    floatGray = 255.0*floatGray;
//    floatGray.convertTo(res, CV_8UC1);
    normalize(floatGray, res, 0, 1, NORM_MINMAX, -1, mask);
    return res;
}


/*
 * This function auto thresh the input image. using binary search.
 * THe final result is an CV_8U image, with write/black ratio
 * approximately at $target.
 */
Mat autoThresh(Mat input, double target, InputArray _mask, int numIt) {
    Mat mask;
    _mask.copyTo(mask);
    CHECK(mask.empty() ||
          (mask.rows == input.rows && mask.cols == input.cols)) <<
                                                                "size of mask and input image must be same";
    CHECK(input.channels() == 1) << "Must be a 1-channel image";
    if (mask.depth() != CV_8U)
        mask = myIm2Int(mask);
    double thresh_h, thresh_l = 0;
    if (input.depth() == CV_8U)
        thresh_h = 255;
    else
        thresh_h = 1;
    double allPics;
    double thresh;
    if (mask.empty()) {
        allPics = input.rows * input.cols;
    } else {
        mask = myIm2Double(mask);
        allPics = sum(mask).val[0];
    }
    for (int i = 0;; i++) {
        thresh = (thresh_h + thresh_l) / 2;
        double s = sum(myIm2Double(input > thresh).mul(mask)).val[0];
        if (s / allPics < target) {
            thresh_h = thresh;
        } else {
            thresh_l = thresh;
        }
        if (i == numIt) {
            break;
        }
    }
    return input > thresh;
}

Mat getCentrelMask(Mat input) {
    CHECK(input.depth() == CV_8U) << "Only CV_8U format supported";
    CHECK(input.channels() == 3) << "Must be three channel image";
//    myImWrite(input,"oriimg", id);
    int sizer = input.rows, sizec = input.cols;
    resize(input, input, Size(), 0.5, 0.5);
    Mat mask = create_mask(input);
    vector<Mat> channels;
    split(input, channels);
    Mat green = channels[0];
    medianBlur(green, green, 21);
//    green = grow_boundary(green,mask,21,50);
//    Mat binaryMask = (green > 0) / 255;
//    green = green.mul(binaryMask);
//    green += (1-binaryMask) * 70;

    green = myIm2Double(green);
//    myImWrite(green,"before norm", id);
    Mat norm;
    norm = green;
//    norm = localNorm(green,90,300,mask);
    double vMax, vMin;
    minMaxLoc(norm, &vMin, &vMax);
//    myImWrite(norm,"after norm",id);
    Mat result = norm > (vMax * 0.6 + vMin * 0.4);
//    Mat result = autoThresh(norm,0.003,mask,10);
//    result = find_roundObject(result, 50);
    CHECK(result.channels() == 1) << "Error return pic channel";
    CHECK(result.depth() == CV_8U) << "Error return depth";
//    myImWrite(result,"center",id);
    Mat resizeToOriginalSize;
    resize(result, resizeToOriginalSize, Size(sizec, sizer));
    return resizeToOriginalSize;
}


/*
 *  This function returen a mask indicating the optic disk and the center
 *  bright area that should not be regarded as area of leakage.
 *
 *  Retuned mask is a 1-chnnel and cv_8u type image.
 */
Mat getOpticDisk(Mat input) {
    CHECK(input.depth() == CV_8U) << "Only CV_8U format supported";
    CHECK(input.channels() == 3) << "Must be three channel image";
    int sizer = input.rows, sizec = input.cols;
    resize(input, input, Size(), 0.2, 0.2);
    Mat mask = create_mask(input);
    vector<Mat> channels;
    split(input, channels);
    Mat green = channels[2];
    medianBlur(green, green, 21);
    green = grow_boundary(green, mask, 21, 50);
    Mat binaryMask = (green > 0) / 255;
    green = green.mul(binaryMask);
    green += (1 - binaryMask) * 70;

    green = myIm2Double(green);
//    myImWrite(green,"before norm", id);
    Mat norm;
    norm = localNorm(green, 90, 300, mask);
//    myImWrite(norm,"after norm",id);
    Mat result = autoThresh(norm, 0.035, mask, 10);
    result = find_roundObject(result, 50);
    CHECK(result.channels() == 1) << "Error return pic channel";
    CHECK(result.depth() == CV_8U) << "Error return depth";
//    myImWrite(result,"opdisk",id);
    Mat resizeToOriginalSize;
    resize(result, resizeToOriginalSize, Size(sizec, sizer));
    return resizeToOriginalSize;
}


typedef unsigned char uchar;

Mat find_roundObject(Mat img, int radius) {
    static int counter = -1;
    counter++;

    Mat label, stats, centroids;
    CHECK(img.channels() == 1);
    CHECK(img.depth() == CV_8U);
    connectedComponentsWithStats(img, label, stats, centroids, 4, CV_32S);
    Mat oimg(img.rows, img.cols, CV_8UC1);
    oimg = 0;
    Mat templateImg(radius * 2, radius * 2, CV_8U, Scalar(0));
    circle(templateImg, Point(radius, radius), radius, Scalar(255), -1, FILLED, 0);
    Mat matchResult;
    matchTemplate(img, templateImg, matchResult, CV_TM_CCORR_NORMED);

    int sumArea = 0, area = 0;
    map<int, float> similarity, pixelCount;
    for (int i = 1; i < stats.size[0]; i++) {
        sumArea += stats.at<int32_t>(i, CC_STAT_AREA);
    }
    area = (int) (sumArea * 0.2);
    for (int i = 1; i < stats.size[0]; i++) {
        if (stats.at<int32_t>(i, CC_STAT_AREA) < area)
            continue;
        // calculate similarity of area to round object
        int xmin = stats.at<int32_t>(i, CC_STAT_TOP);
        int ymin = stats.at<int32_t>(i, CC_STAT_LEFT);
        int xmax = stats.at<int32_t>(i, CC_STAT_TOP) + stats.at<int32_t>(i, CC_STAT_HEIGHT);
        int ymax = stats.at<int32_t>(i, CC_STAT_LEFT) + stats.at<int32_t>(i, CC_STAT_WIDTH);
        pixelCount[i] = 1;
        similarity[i] = (float) matchShapes(templateImg, img.rowRange(xmin, xmax).colRange(ymin, ymax),
                                            CV_CONTOURS_MATCH_I1, 0);
    }

    int selectedArea = -1;
    float maxSimiliarVal = 100;
    for (auto i : similarity) {
        i.second /= pixelCount[i.first];
        if (i.second < maxSimiliarVal) {
            maxSimiliarVal = i.second;
            selectedArea = i.first;
        }
    }
    for (int x = stats.at<int32_t>(selectedArea, CC_STAT_TOP);
         x < stats.at<int32_t>(selectedArea, CC_STAT_TOP) + stats.at<int32_t>(selectedArea, CC_STAT_HEIGHT);
         x++) {
        for (int y = stats.at<int32_t>(selectedArea, CC_STAT_LEFT);
             y < stats.at<int32_t>(selectedArea, CC_STAT_LEFT) + stats.at<int32_t>(selectedArea, CC_STAT_WIDTH);
             y++) {
            if (label.at<uint32_t>(x, y) == selectedArea) {
                oimg.at<uchar>(x, y) = 255;
            }
        }
    }
    return oimg;
}


